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Abstract. We study stationary axisymmetric configura- 
tions of a star model consisting of two barotropic flu- 
ids, which are uniformly rotating at two different rotation 
rates. Analytic approximate solutions in the limit of slow 
rotation are obtained with the classical method of Chan- 
drasekhar, which consists of an expansion of the solution 
in terms of the rotation rate, and which is generalized to 
the case of two fluids in order to apply it to the present 
problem. This work has a direct application to neutron 
star models, in which the neutron superfluid can rotate 
at a different speed than the fluid of charged components. 
Two cases are considered, the case of two non-interacting 
fluids, and the case of an interaction of a special type, cor- 
responding to the vortices of the neutron superfluid being 
completely pinned to the second fluid. The special case 
of the equation of state P ex p^ is solved explicitly as an 
illustration of the foregoing results. 
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1. Introduction 

More than 30 years after the discovery of the pulsar phe- 
nomenon and its identiflcation with rotating neutron stars 



(Gold 1968), there exists today a considerable body of ob- 
servational data (Lyne & Graham-Smith 199S), but also 



still a number of uncertainties and open questions about 
the theoretical model for pulsars, mainly due to the ex- 
tremely dense (and theref ore p oorly known) state of mat- 
ter implied (Glendenning 1997 ). 

One of the generally agreed characteristics of neutron 
stars is the existence of a superfluid neutron component. 
This is not only predicted by calculations from nuclear 



physics (Ainsworth et al. 1989), but also agrees with ob- 
served features of pulsar behavior, like the very long re- 
laxation times, from months up to years, after a glitch 
(a sudden increase of angular velocity of the order of 
Afl/il < 10^^). All the charged parts of a neutron star 
(nuclei, protons and electrons) can be treated as a single 
normal fluid, and are predicted to be "locked" together 



in a state of corotation (Easson 1979) on sufficiently long 
timescales. In contrast, the neutron superfluid can have a 
different rotation even on very long timescales, and so one 
is naturally led to consider a neutron star model consist- 
ing of two independent fluids, an approach that was first 



adopted by Baym et al. (1969). This model has since been 



the basis of our understanding of the glitch behavior and 
the subsequent post-glitch relaxation observed in pulsars 



(Anderson & Itoh 1975; Alpar et al. 1984; Sedrakian et al. 
p^95b[ 



Link & Epstein 1996) 



Apart from being inviscid and therefore forming an in- 
dependent fluid component, a superfluid is moreover con- 
strained to be in a state of irrotational flow, and conse- 
quently its rotation can only be achieved by the presence 
of quantized vortices. These vortices will interact with the 
fluid of charged compo nents (Feibe lman |l97l| ; Sauls et al. 
1982| ; Epst ein & Baym |l988| ; Jone s |l99C| , |l99l[ ; Link & Ep- 



stein 1991 ; Sedrakian & Sedrakian 1995a ), giving rise to an 
effective friction force on a moving vortex, and they can 
even be completely pinned to the Coulomb lattice of nuclei 
that forms the crust of the neutron star. A consequence is, 
that the vortices will not corotate with the superfluid and 
will therefore be subject to the Magnus force orthogonal 
to their relative velocity with respect to the superfluid. 
These forces will balance each other, which leads to an 
effective interaction between the two fluids. 

The long-term slowdown of the neutron star's rotation 
rate, which is caused by the loss of energy in form of elec- 
tromagnetic radiation, has many important consequences. 
The global slowdown tends to decrease the ellipticity of 
the equilibrium shape of the neutron star. This leads to 
the buildup of stress forces in the solid crust, which can get 
suddenly released in form of a starquake. This has been 



proposed by Ruderman (1969) as one of the first models 
in order to explain glitches, and has since been a subject 
of great interest, directly as a model for glitches (Baym 
& Pines 1971 Heintzmann et al. 1973; Ruderman 1991 



Link et al. 1998), or at least as a trigger for some other 



glitch-mechanism via the energy liberated in such a star- 



quake event (Link & Epstein 1996). Another aspect of the 



glob al slo wing down has been pointed out by Reiseneg- 
ger ( 1995| ): the decrease of the centrifugal force leads to a 
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global compression of the neutron star matter (consisting 
of neutrons, protons and electrons). But the equilibrium 
composition (with respect to /3 reactions) of this plasma 
depends on the density, and so a global compression drives 
the plasma out of equilibrium. This has some possibly ob- 
servable consequences, e.g., on the emission of neutrinos 
and on the evolution of the temperature of neutron stars. 
These consequences have been examined from the 
point of view of a global slowdown of the whole neutron 
star, but it has to be noted that in the two-fluid model, 
it is primarily the fluid of charged components that gets 
slowed down, while the superfluid neutrons will signifi- 
cantly lag behind and continue to turn at a faster rota- 
tion rate. It has been remarked recently (Carter et al. 



199£), that this could lead to a new mechanism to induce 



stress forces in the crust, due to an increasing deficit of 
centrifugal buoyancy. The model for the driven deviation 
from chemical equilibrium also has to be refined according 
to the two-fluid picture. Not only is there a global com- 
pression, but also a relative displacement of the two mass 
distributions with respect to each other, as the difference 
of their rotation rates increases. For example, when the 
two fluids have been in /3 equilibrium in the state of coro- 
tation, the slowdown of one fluid changes its ellipticity and 
therefore moves volume elements of that fluid to regions 
with a different equilibrium composition, so that they are 
no longer in a state of equilibrium with the second fluid. 
The purpose of this paper is to study the consequences 
of the two fluids having different rotation rates on the 
mass distribution of the star. Even in the case of a single 
rotating, self-gravitating fluid, it is impossible to obtain 
exact analytic solutions, and one has to rely either on nu- 
merical treatments or on analytic approximations (e.g., 
see Tassoul 19781). In the present work we will develop 



a generalization of the analytic approximation of Chan- 
drasekhar (Chandrasekhar [l935 ; Tassoul 1978) to the case 
of a barotropic two-fluid star. This method consists of an 
expansion of the rotating solution around the static solu- 
tion in terms of the rotation rate. Using this method, we 
will obtain an expression for the stationary mass distri- 
bution of a barotropic two-fluid star up to second order 
in the two rotation rates. The obvious limitations of this 
approach are that the rotation rates have to be small com- 
pared to their "natural" scale, and that both have to be 
of the same order of magnitude. These conditions are in 
general satisfied in the case of neutron stars. The fact that 
we considered stationary solutions is no real restriction ei- 
ther, as the slowdown of pulsars takes place on very long 
timescales. Therefore it should be possible to describe it 
as a quasi-stationary process, passing through a series of 
stationary states. 

The plan of this paper is the following. In section pi 
we define the Newtonian general model of a barotropic 
two-fluid star, and in section || we further specialize this 
general model in the context of neutron stars. In section 
y we generalize and apply the classical method of Chan- 



drasekhar to this two-fluid star, which allows us to reduce 
the problem to a set of ordinary differential equations. Sec- 
tion H is devoted to the boundary conditions necessary to 
obtain the complete solution, which is given in section pi 
Section ^ is concerned with some consequences of the so- 
lution, like the change in ellipticity and moment of inertia. 
In section H we discuss an effect that we call "rotational 
coupling" , which is the fact that changes of the rotation 
speed of one fluid influence the rotation of the other fluid 
via the gravitational potential, even if the two fluids are 
supposed to be strictly non-interacting. Section ^ gives 
an illustration of the foregoing results in the completely 
analytically solvable case of a special polytropic equation 
of state. Section |lC| summarizes this work. 

2. The two fluid model 

We want to describe a star consisting of two independent 
fluids in Newtonian gravitation. We distinguish a fluid de- 
noted by the subscript c, that will represent the globally 
neutral fluid of charged components of a neutron star (nu- 
clei of the crust, protons and electrons), and a fluid de- 
noted by the subscript s, that will describe the superfluid 
of free neutrons. We will also refer to the fluid of charged 
components as the "normal fluid" , as opposed to the su- 
perfluid. So the basic description of our model consists of 
the Eulcr equations for the two fluids: 

(1) 

Pc {dtvl + vi^jvl) = -V^Pc - PcVV + fi , 

where dt denotes the partial derivative with respect to 
time, Pa, Pa, Va ^'^d fa ^^'^ the respective mass density, 
pressure, velocity and force per volume of each of the two 
fluids, and a is the "chemical index" (a = s,c). is the 
gravitational potential, which is related to the total den- 
sity p = pc + Ps by Poisson's equation 



V^ == AnGp, 



(2) 



where G is Newton's constant. 

We consider only stationary, axisymmetric configura- 
tions, with the two fluids rotating uniformly with respec- 
tive angular velocities fij. and i^s, i.e., Va = J^a x f. In the 
subsequent analysis we work with dimensionless quanti- 
ties, measuring length scales in units of the radius R, den- 
sities in units of the central density po of the non-rotating 
conflguration and time in units of l/y/AirGpo . Table [^ 
shows a summary of the employed fundamental and de- 
rived units. In order to avoid unnecessary complications 
of notation, we will in the following keep the same sym- 
bols for the dimensionless variables, with the exception of 
the rotation rates Hq,, which we will now denote £„. This 
is in order to emphasize the fact that we are considering 
slow rotations with respect to the natural scale of fl (see 
Table |l|; this scale is in general still bigger than the Ke- 
plerian rotation rate f2^ = 47rGp/3, where p is the mean 
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Table 1. The system of the chosen "natural" units, R and 
Po are respectively the radius and the central density of 



the non-rotating configuration 
Quantity 



Unit 



Length 


R 


Density 


Po 


Time 


l/V47rGpo 


Frequency 


^AttGpo 


Mass 


PoR' 


Moment of Inertia 


poR' 


Gravitational Potential 


A-wGpoR^ 


Pressure 


AtiGpIR" 


Angular Momentum 


V^TvGpf/'R" 


Force/Volume 


inGp'oR 



density), and therefore £„ represents a small parameter, 
i.e., Sa <. 1. 

The rescaled Euler equations M) take the form 



-V^P. + V*( 






Pc 2 



Ps 

1 

Pc 



(3) 



-fl 



where w is the cylindrical radius, and Poisson's equation 
(0) in the new variables reads 



V^(/. = p. 



(4) 



The fundamental assumption in our treatment is that 
each of the two fluids obeys a barotropic equation of state 
(EOS), i.e.. Pa — Pa{Pa)- This allows us to write the 
terms V^Pa/pa in ® as the gradient of a function —^a, 
say, that is defined as 



Ipa 



' dp 

Pa{p) 



(5) 



As we will work in the approximation of T = 0, the quan- 
tity —tpa is equivalent to the enthalpy per mass unit and 
to the chemical potential per mass unit, and in our subse- 
quent analysis it will play the role of an effective potential. 

3. The two— fluid model for neutron stars 

In the present section we will specialize the general two- 
fluid model of the previous section to the case of a neutron 
star. The "normal" fluid of charged components (c) is sup- 
posed to be corotating with the crust on short timescales, 
because of the strong magnetic field that "locks" all 



charged components to the crust (Easson 1979). The in- 
dependent component (s) is a perfect fluid that coexists 
with the normal fluid without any viscous interaction, but 
we will allow for an indirect interaction via the super- 
fluid vortices. We neglect all magnetic and thermal influ- 
ences, as we are mainly interested in the effects of rota- 



tion. While the assumption of uniform rotation is proba- 
bly quite realistic for the normal fluid, the superfluid neu- 
trons could perfectly well be in state of differential rotation 
{W^fls ^ 0), even under the condition of stationarity, but 
for simplicity we will assume it to be uniformly rotating. 

As we are interested in stationary solutions, we will 
also neglect the external forces acting on the neutron star, 
which, for isolated neutron stars, are due to electromag- 
netic radiation and lead to the long-term slowdown of the 
rotation rate of the crust. This approximation is easily jus- 
tified, as the timescales of mechanical displacements of the 
neutron star matter due to a change in rotation is much 
shorter than the typical slowdown-timescale fi/fi, which 
is of the order of 10^ years. 

Part of the normal fluid, namely the solid crust, is not 
really a fluid, but we could still approximately describe 
it as a fluid subject to anisotropic volume-forces, namely 
the stress forces due to the solidity. This means that the 
force density /* acting on the normal fluid would not only 
consist of the opposite interaction term —fl, but also of 
an extra term fl due to the anisotropic stress forces, i.e., 
we would have 



ri ri ri 

J c J a J s 



(6) 



The fact that there is no temperature-dependence in the 
bulk EOS is an excellent approximation in the neutron- 
star context, as the actual temperatures (for not extremely 
young neutron stars) are some orders of magnitudes be- 
low the Fermi temperature. Additionally, as we assume 
two independently conserved barotropes, we also neglect 
possible "chemical interactions" between the two fluids 
via P reactions, which transform neutrons into protons 
and electrons and vice versa {n ^^ p + e + De) . But the 
nature of the involved f3 reactions in neutron star matter 
(namely, indirect Urea) seems to be rather slow, i.e., the 
chemical equilibration timescales are of the order o f sev- 
eral years for not very young neutron stars (Haensel 1992 ) 
and therefore much longer than the dynamical timescales 
under consideration. So the above approximation should 
be rather viable as long as we do not consider evolutions 
on very long timescales, where inevitable effects of trans- 
fusion would have to be included in the analysis (e.g., see 
Langlois et al 



1998) 



We still need to specify the nature of the interaction 
force /*. The conditions of stationarity and different rota- 
tion rates do not allow a dissipative interaction between 
the superfluid vortices and the normal fluid, so we are 
basically left with two possible types of interaction, the 
case of completely pinned vortices, e.g., as obtained by Ep- 
stein and Baym ( 1988| ), and the case of quasi-/ree vortices, 
as suggested by the results of Jones (1991). The pinned 
case should still be a good approximation even if vortex- 
creep is effective (that is, the vortices jump from pin- 
ning site to pinning site, as they are pushed by the Mag- 
nus force), whenever the creep-timescale is long compared 
to the dynamical timescale, so that the quasi-stationary 
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mass distribution in the creep case should not differ from 
the pinned case. The pinned case leads to an interaction 
caused to the Magnus force acting on the vortices, which 
is given by 



fli = Ps{es~Sc)es'^'tu'^ ■ 



(7) 



This supposes a parallel lattice of vortices. We will follow 
this common assumption, which has been shown to be 
valid under certain conditions by Ruderman and Suther- 



land (1974). In the free case we have 

/: = o, (8) 

so we can treat the two cases (0) and (^) together, writing 

/: = -^p/m, (9) 

where the "pinning switch" Sp is 1 in the pinned case and 
in the free case. 

We arrive at the following form for the two Euler equa- 
tions ra): 



'^ps+(f)- 



1 



e„n7 



SpEsiss - ec)vj = 



(10) 



applied Magnus force. On the other hand, contrary to a 
fluid it is able to keep the vortices from moving relative 
to the crust. 

The condition of constant «;(r) = k does not restrict 
the choice of the EOS of both fluids, but only fixes the 
EOS of the second fluid with respect to the chosen EOS 
for the first fiuid by the relation 



Pc{Pc) = -Ps{nPc 
K 



(13) 



In the following we set fl = and postpone the difficult 
problem of including anisotropic stress forces to future 
work, so we restrict our analysis to the two above men- 
tioned completely "fluid" cases: 



(i) free vortices {Sp = 0) 
(ii) pinned vortices {5p = 1) 

(with the EOS subject to (|lj 
is a constant) 



such that ps / Pc 



From equations 
potentials 



(M) and fl = we obtain the effective 



1 



Pc 



where we have defined 
Pc '{r) 



(11) 



(12) 



Pa being the zeroth order density distributions, that is, 
of the non-rotating configuration. 

We see that the right-hand side of (|ll|) has to be the 
gradient of some scalar function. Looking at the pinning 
term (containing Sp) of this equation, we see that this 
term alone can in general not be written as a gradient, 
because of the factor K{r). This shows that in general the 
pinning force cannot be compensated without the presence 
of the anisotropic stress force /^ , which is provided by the 
solidity of the crust, as has already been noticed in the 



literature (e.g., Ruderman [l99lD . 

There is however a special case that has the advantage 
of being analytically tractable, where the pinning force 
can be compensated by the gradient force on the left- 
hand side alone, without including any stress forces. This 
is obviously the case when K(r) is a constant. As with 
our preceding assumption of uniformity of fig , this case is 
not necessarily realistic for neutron stars, but it is still of 
interest since it provides qualitative insight in the behavior 
of the system in the pinned case. It corresponds to the 
limiting case of a very ductile crust that does not develop 
any notable shear stress and deforms like a fluid under the 



V-s = - — (e' - 2Spe,ies - £c)) + C, , (14) 

2 

ipc = '/' - ^ (Ec + ^nSpEsies - £c)) + Cc , (15) 

where the Ca are constants in space, but they may depend 
on the rotation rates Eq . One can see that the pinned case 
{Sp = 1) introduces mixed terms Ss£c, while in the free 
case {Sp = 0) the only non-zero terms are the diagonal 
ones, that is e^. 

The pressure Pa should be a monotonic function of 
density pa, and so we see from (^ that V'a should also be 
a monotonic function of pa ■ This relation can therefore be 
globally inverted, so that the density pa can be uniquely 
written as a function of the effective potential ipa in the 
form 



Pa = Pa{i^a) 



(16) 



a relation that will be important for the subsequent anal- 
ysis. 

4. Generalized Chandrasekhar expansion 

It will be convenient, in order to obtain more compact ex- 
pressions, to introduce a matrix notation in the fluid in- 
dices. One will effectively recover the usual Chandrasekhar 
type of terms know from the case of one fluid (see Chan- 
drasekhar 1933 ;Tassoul 197§|) , with the scalar perturba- 
tion quantities replaced by symmetric 2x2 matrices. We 
write the effective potentials as follows: 



Ipa = (/>■ 



— -e • Za ■ £ + Ca (s) , 



(17) 
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where the "centrifugal" matrices Z^ are defined as 



^■^ = ' \)-'5M -1 0^ 



Zc ^ \ 1 ) + '«'5p ( _i 



(18) 



By writing M we indicate that the quantity M is a sym- 
metric 2x2 matrix in the fluid indices with components 
M"^, and e is the vector with components Eq. 

Following the standard method of Chandrasekhar, we 
expand all quantities up to second order in the rotation 
parameter e around the non-rotating configuration. Be- 
cause of the symmetry under parity, i.e., e — )■ — e, there 
can be no terms of first order in e. The second-order term 
is a quadratic form in e and therefore the definition of the 
coefficient matrix is ambiguous. We can fix this ambiguity 
by the additional condition that the matrices occurring in 
the expansions have to be symmetric. 

We work in spherical coordinates r and u = cos(6') 
(where, of course, 9 is defined with respect to the axis of 
rotation) and so for the fluid densities pa {r, u) this expan- 
sion reads 



Pa {r, u) = p^"^ (r) + 5pa (r, u) with 



e-pa-e. (19) 



We expand the other quantities 0, Cq and p in the same 
way, with the respective second order coefficient matrices 
0, Ca and p (where of course p = ps + Pc)- 

It is important to note that the additive constants Ca 
depend in general on the rotation rates e. We can absorb 
the additive constant Ca into the definition of Pa{'4'a), 
so for convenience we can set Ca — 0, but we have to 
keep track of the 0{e^) correction e • Cq • e. 

In order to obtain the relations between pa and to 
second order in e, we expand pa{ipa) around the non- 
rotating configuration i/jq ' — 0^°^ : 

pMc) = Pa(<^^"^)-fca e- U - ^Za + Co) •e+0(e^) , (20) 
where 

_ dpailp) 



ka 



dil; 



(21) 



,^(0) 



Order by order comparison between ( pO| ) and (M) together 
with the condition of symmetric matrices leads to the iden- 
tifications 



Pa ~^ ^f^a 



~^a + ^a 



(22) 



which further allows us to write the "structure function" 
ka simply as 



kair) = - 



^ (0) 
dpa 



(23) 



p=-kd+^K{r)-D{r), 

where we have defined k = kr + kg and the matrices 



The total density perturbation coefficient is found from 
l2|) to be 



(24) 



K{r) = - iksZg + keZA and D{r) = kgC^+k^Cc .(25) 

Surprisingly, the matrix K is found (using the definitions 
of K and ka, (O) and (p3)) to be the same in the free (i) 
and the pinned (ii) case, namely. 



3 I fee 



(26) 



Inserting ( p4|) into Poisson's equation (|j), one finally ob- 
tains the partial differential equation for the second order 
corrections (p of the gravitational potential. 



V^^ + /c = -^ K{r) - D{r) 



(27) 



Using the decomposition of (/)(r, u) in the orthogonal basis 
of Legendre polynomials, we can reduce this partial differ- 
ential equation to an infinite series of ordinary differential 
equations. We write 



^w) =^ P2i{u)(j)2i{:r 



(28) 



i=0 



where we only need to sum over Legendre polynomials 
with even index, assuming equatorial symmetry. Using the 
well known differential equation for the Legendre Polyno- 
mials, the Laplace operator acting on (j) is seen to reduce 
to 



V20(r,7/)-^P2i(wp2i 02^(0, 



1=0 



where the differential operator I?„ is defined as 



2?„ = 



_ ^2 2d n{n + 1) 



dr^ r dr 



(29) 



(30) 



Using the orthogonality property of the Legendre polyno- 
mials together with the fact that 3vu'^ = 2r'^{l — P2iu)) 
leads to the following series of ordinary differential equa- 
tions 



Vo4>o + k4>o = +r^K -D, 
I?202+fc02 = -r^K, 
V2i^2i + k^2i =0 for l>2. 



(31) 



In order to solve these equations, one must specify the 
appropriate boundary conditions, which we consider in the 
next section. 
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5. Boundary conditions 



The first restriction on the solutions of (31) comes from 
the requirement that the 02; should be regular functions 
in r = 0, and therefore the left hand side of the differential 
equation has to be regular in the origin too. This leads to 
the conditions 



E>2/ 



(0) 



02^(0) 



for l> 1, 



for l>0. 



(32) 



The prime stands for derivatives with respect to r. An- 
other boundary condition is obtained by matching the so- 
lution for the gravitational potential inside the star to the 
solution (j)E outside the star. The external solution is nor- 
malized conventionally by lirur^oo (f'E = 0, and satisfies 
^'^4'E = 0. Its expansion in terms of Legendre polynomi- 
als, and up to second order in e has therefore the following 
form: 



..(") 



<i>E(r,u)='—^e_-\y^^P^Au)\ ■e_^0(e^). (33) 

Taking into account the deviation of the star from spheric- 
ity, the surface can be expressed as 



(34) 



i?(ii) = l + e- y^RiiP2i{u)\ •£ + 0(e^), 

where the radius of the non-rotating configuration i?'"' — 
1 in our units (see Table |l|) . The matching conditions are 
given by 



(f>{R{u),u) = (f)E{R{u),u) , 



c^'{R{u),u) = (f>'j^{R{u),u). 



(35) 



The deviation of the derivative normal to the surface from 
a simple radial derivative is of order ©(e"*), so we can 
neglect it. 

Expanding these matching conditions up to second or- 
der and using the fact that 0(")"(1) + 20("''(1) = yields 
the following boundary condition for the (/)2/: 



^^,(l) + (2; + l)02;(l) = O. 



(36) 



It is interesting to note that this condition was found with- 
out ever specifying the actual surface of matching. The R21 
were in fact completely arbitrary apart from the restric- 
tion to be small compared to e^^, such that the devel- 
opment (34) makes sense. This shows that the obtained 
boundary relation for the (f)2i is a rather robust conse- 
quence of the matching to the vacuum solution itself. One 
could in fact find the R21 which specify the actual surface 
of the star up to second order in s in terms of the (f)2i by 
the obvious definition 



p{R{u),u) =0 



(37) 



which then leads to the expression for the surface up to 
O(e^) in the form 



R(u) = 1 ; £• I Vp2/(1)-P2i(^t) I ■£■ 

p(o)'(l)- [^/''^ ^ ''^ ') - 



(38) 



For the individual fluids we can find the Pa — surfaces 
in the same way: 



Rju)^R^°'^~ 



1 



Pa (rCa ) \l=0 



e-ij2p-MR^°^)P2du))-e.m 



It has already been recognized by various authors that this 
type of expansion eventually becomes singular in the vicin- 
ity of the star's surface (see Smith 1975; Tassoul 197§| , and 



references therein). The zeroth order term oi p^'^^ {r)+e-p-e 
obviously becomes zero on the non-rotating star's radius 
r = 1, and so the 0{e^) correction can no longer be con- 
sidered as being small with respect to the zeroth order 
term. Due to this fact the value for p(r, u) is locally valid 
only as long as one stays away from the surface, and so the 
definition of R21 via ( |37| ) seems rather unreliable. There- 
fore it is important that the boundary condition (36) does 



not depend on the actual form of the boundary surface. 

We note that for the case ^ = we still need two more 
conditions in order to fix all the 4 free parameters of the 
solutions 4'Q{r) and pa^o (C'q and the two free parameters 
for a solution of a differential equation of second order). 
These conditions are obtained by invoking the requirement 
of mass conservation for each of the two fiuids: 



d^xpa,{r,u)= / d^xp^°\r). 



(40) 



.(0) 



The fact that pa (r) vanishes in Ra leads to 



(0) 



d^xpa{r) 



-he 



d'xp^^\r) 



d xpc 



+ Oie'). (41) 



Because of the orthogonality property of the Legen- 
dre polynomials and Po(''^) = 1; f^ny integral of the type 
j_^duP2i{u) vanishes for / 7^ 0, so that the condition of 
mass conservation simply reduces to 



drr'^Pa.oir) = 0. 



(42) 



As mentioned by Heintzmann et al. (1973) in the case of 
one fluid, the integral constraint of totalmass conservation 
can be reduced, with the help of Poisson's equation (0), 
to a differential boundary condition on (J)q, namely. 



?o(l) = 0, 



(43) 



but in the case of two fluids considered here, we still have 
to use one of the two integral constraints ([42|), in order to 
fix the second constant. 
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If one wanted to consider a transfusive type of model 



(see Langlois et al. 1998), one would effectively have only 
( p^ ) and would still need some other prescription in order 
to fix the remaining constant, and thereby the respective 
transfusive mass transfer between the two fluids. 



poles P2{u) = +1, so we get the general expression for the 
ellipticity: 



(-s-PaARa) -e) 






JO)' 



(i?i"^) 



(45) 



6. Formal Solution 

The prescription of the boundary conditions not only com- 
pletely specifies the solutions of our series of ordinary dif- 
ferential equations (pi]), but it even restricts nearly all of 
them to be zero. For / > 2, (j)2i is given by a homoge- 
neous differential equation of second order, subject to the 
boundary conditions ( p2| ) and (36). Only one of the two 
fundamental solutions can be chosen to be regular in the 
origin, so we have the freedom of only one multiplicative 
constant in order to satisfy (pq) , which can in general only 
be zero. All the solutions arc trivial whenever the differ- 
ential equation is homogeneous. This is the case for all the 
02/ with I > 2, but also for those matrix-elements in the 
cases I = and ^ = 1, for which the inhomogeneous term, 
that is the corresponding matrix-element of K and I?, is 
zero. 

Taking a look at the elements of K in the free (i) and 
the pinned (ii) case (Eq) , we immediately see that the only 
non-trivial solutions will be the diagonal ones, 0"" and 
(p'^'^. So the formal solution of the problem consists of the 
following density perturbation coefficients (see (|2S 



(44) 



Pafii"^) = -ka{r) [ (poir) - —Za + Ca 



Pa,2{l^) = ~ka{r) ( 02 (r) + —Za 



with the (f)o{r) and 02(?') solutions of (|3l|), subject to the 
conditions of regularity (|32| ) , continuity with the external 
potential (pq), and mass conservation ([4^). 



7. Expansion, Ellipticity and Moment of Inertia 

We will now discuss some of the consequences of the ob- 
tained formal solution up to second order in e for the 
densities pa (r, u) . For simplicity, we restrict our attention 
in this section to mass distributions Pa with a simply con- 
nected topology, that is to say, which possess only one 
boundary surface for each fluid, namely the outer surface, 

and so we have pa (Ra ) < 0. From the expression for 
the respective boundary surfaces ( p9| ) we see that there 
is a uniform expansion of the fluid as a whole of amount 

e ■ i^Pa.oiRa )l Pa {Ra )) ' £, and superposed on this a 
term proportional to P2{u), which leads to the ellipticity 
of the surface. At the equator P2{u) = —1/2 and at the 



From ([44| ) and the regularity condition 02(0) = 0, we see 
that p2(0) — 0, and so the relative change of the central 
density is given by 



PaiO) e-Pa,o(0)-e 



pi^O) 



p'^\0) 



(46) 



We should also expect the volume to change, by an amount 
given by 



y„ = yW+4^(i?i"))^ 



£■ Pafi{Ra ) ■ £ 



Pa (rCa , 



(47) 



Finally, we write the change of the two moments of inertia 
in the form 



where /„ is given by 

d XTU Pa{r, u) . 



(48) 



(49) 



We note that the integral is done only over the unper- 
turbed, spherical volume Va , because the corrections due 
to the form of the boundary surface are of order 0{e^), 
which is due to the same cancellation as has already been 
encountered in the density integration (|4l|). Further eval- 
uation leads to 



''-il -' 



drr 



Pafiir) - gPa,2(0 



(50) 



8. Rotational coupling 



In this section we investigate a consequence of the depen- 
dence of the moments of inertia on the rotation rates e, 
which is expressed in equation (Hq) . The moment of iner- 
tia of one fluid also depends on the rotation of the second 
fluid, which leads to what can be called "rotational cou- 
pling" . This effect is still present in the free case, where 
the only way the two fluids communicate with each other 
is via the gravitational potential 0: changing the rotation 
rate of the fluid a changes its mass distribution pa and 
therefore also 0, which in its turn will change the mass 
distribution of the second fluid pp. As we saw above, this 
effect takes place on the order 0{e'^). 

Let us consider the angular momentum, which in our 
units (see Table 1^) is given by 



La^ (4°^ +e-T^-e) £^ +Oie^) 



(51) 
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If we want to express the rotation rates Eq. in terms of the 
angular momenta Lq,, it suffices to invert this relation and 
we obtain 



= £(1) 



.(1) 



f- 



I. 



r(l) 



r(0) 



(52) 



where we have defined the first order rotation rate by 
(1) ^ L^ 



r(0) 



(53) 



which is the rotation rate for a given angular momentum 
-Lq, if we kept the mass distribution fixed to the value of 



the non-rotating case. We see that in (52), at the order 



0{e^), we were allowed to replace £« by e^ . This is the ex- 



plicit relation for £„(£« ^ Ec^Oj or equivalently £a{Ls, Lc). 



-(1) ^(l)^ 

Here we see again the effect of the rotational coupling be- 
tween the two fluids, namely, the change of the rotation 
rate of one fluid if we change the angular momentum of 
the other fluid. This mutual dependence explicitly reads 
as 

^"" ''^" "' ^ ^(i)^<5„.-2£Wf/„.£Wl .(54) 



dL 



P 



{iL' 



= ( /1°) 



.(1) 



3q/3- 



i"( 
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The effect is of order 0{e'^) and its actual importance is 
determined by the coefficient la /la , which depends on 
the EOS. 

Let us take a look at a particular case, where we change 
the angular momentum Lc without changing L^, corre- 
sponding to what happens in a real neutron star, for ex- 
ample when we consider the loss of angular momentum of 
the normal fluid due to electromagnetic radiation. In this 
case we can express the change of angular velocity of the 
superfluid with respect to the change of the normal fluid 
as 



dsc dsc/dLc 



= -2eli) 



r(0)- 



.(1) 



0{e'). 



(55) 



9. Exact solution 

for the polytrope P oc p^ 

In the previous sections we have obtained formal solutions, 
and all quantities have been expressed in terms of 0o a-nd 
02 , which satisfy the differential equations (|3l| ) . The pur- 
pose of this section is to consider a special case for which 
these equations can be explicitly solved, and that is the 
case of the two fluids obeying a polytropic EOS of the 
type 



Pa^—- 



(56) 



where for the moment the ka is just a fluid-specific posi- 



tive constant. We can see that the two EOS (56) satisfy the 
relation (O) with k, = ks/kc, so we can study the free and 
the pinned case for this special EOS. The solutions pa (r) 



„(o) 



for the non-rotating case will satisfy pi — upc ■ This 
relation tells us that both fluids share the same bound- 
ary surface, which is therefore the star's surface, and so 

We start by the zeroth order approximation, that is 
the non-rotating configuration of the two-fluid star. The 
equation of hydrostatic equilibrium in the non-rotating 
case reads as 



1 

P«) 



V'Pa = -V 



- _Y7^^(0) 



and for the EOS ( p6| ) it has the solution 



(57) 



(58) 



Using the definition (|23| ) of the structure function ka(r) 
we see it is equal the constant ka defined in the EOS ( p6[) 
For the total density we find 



P 



(0) 



^fc0(") + C(°) . 



(59) 



Putting this expression into Poisson's equation (0), we 
recover the same Lane-Emden equation we would get for 
one polytrope of the form P = (? jlk^ namely, 



V2p(o)(r) + fcp('')(r)=0, 



(60) 



even if the combined system of the two fluids can not be 
described as a barotrope at all, i.e., P{ps,pc) = Ps + Pc 
can not be written as a function of p alone. The above 
equation, subject to the boundary condition p^^'{0) = 1 
(which is due to our choice of units), has the following 
solution: 



p^">{r) 



SlU' 



(rVk) 



(rVk) ' 
which implies that 



for r < 1 , 



(61) 



(62) 



This is not too surprising, as it is well known that in the 
case of a static polytrope with polytropic index 2 there 
exists a simple proportionality relation between the star's 
radius R and the coefficient A:, the radius being in fact de- 
generate with respect to the star's mass. As we are working 
in units where i? = 1, this also fixes the numerical value 
of k. Due to the proportionality relation ps = upc and 
K = ks/kc, we obtain for the respective densities 



(0) ka sin(r7r) 



r-K 



(63) 



We come now to the corrections of order ©(e^), deter- 
mined by the coefficients 0o and 02 that are the solu- 
tions of (|3l|). The regular homogeneous solution is found 
in terms of the spherical Bessel function j„ [x] to be 



<?2i {r) = AJ2i (r7r) . 



(64) 
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Particular solutions are found by inspection, and so we 
obtain the exact solution to (5l^) in the form 



r.^ 7 ■ r ^ ^ f ^ 6\ D 

k{r) = Aojoirn) + ^ k - — ]-—, 



'- - K 

02 (r) = A2J2{rn) -r'^ , 

Mr) =0, for l>2, 



(65) 



where the remaining constants Aq, A2 and D are to be de- 
termined by the boundary conditions (Q) and (p^), which 
finally yields 



K 



Mr) = — (2io(r7r)+r2-l) 



K 



'2W = -^ (5j2(r-7r) -' 



(66) 

(67) 



Inserting the obtained cj) into the equation p^ for the (M) 
and invoking the mass conservation condition ( ^2| ) for the 
individual fluids determines the remaining constants C^. 
For the sake of completeness we will write the complete 
solution (Kh after putting all the pieces together: 



^ / ^ , ] K f^ , ^ 2 3 6 
Pa.o{r) ^ -'^a<^[2jo[rTT)+r 

TT^ V 5 TT 




(68) 



i^ 



PaAr) = -fca -^ — (5 J2(r-7r) - r^) + — -Zq 



while he total density perturbation coefficients p2i can be 
written more compactly, 



Po = -K ( 2 join:) - ^ ) : 



P2 = -bKJ2{rTT). 



(69) 



Using this explicit solution we can evaluate the coefficients 
that determine the rotational coupling (|54| ) discussed in 
section |8[ The integration ( M ) over the explicit solutions 
yields 



(70) 



-j^ = aK + bZa, 



with the coefficients 
9 



3- — 



b = 



7r2 -6 V 5 175; ' 

3^6 



The expression (pq), which applies to the particular case 
where dLg — 0, dLc ^ 0, can now be obtained explicitly 
as 



d£r_ 



-2aKe^^h^P 



in the free case (i), and 
in the pinned case (ii). 



(71) 



(72) 



175(7r2 - 6) 



10. Conclusions 

We have considered stationary, axisymmetric configura- 
tions of two fluids rotating uniformly with different ro- 
tation rates. The analytical method of Chandrasekhar, 
known from the classical problem of a single rotating fluid, 
has been generalized to the two-fluid case. By applying 
this method we have obtained the formal solution of the 
respective equilibrium mass distributions for the two fluids 
in terms of the two functions 00 {r) and 02 (r) , which are 
the solutions of the ordinary differential equations (pi]). 
In order to fully determine these solutions, one needs to 
specify an EOS for the two fluids. The case of the special 
polytropic EOS P ex p^ is solved as an example in section 
0. A genuine effect of the two-fluid model is pointed out 
in section g|, namely, the fact that the gravitational po- 
tential communicates changes in rotation speed and mass 
distribution between the two fluids. 

Further effort would be necessary in order to include 
the effects of solidity of the crust, so that one could an- 
alyze the buildup of stress forces in the crust, including 
the case of pinned vortices in its generality, without the 
present restriction of (|3|) . Further investigations will also 
be concerned with the implications of the present results 
on the deviation from chemical equilibrium and thus heat- 
ing and neutrino emission. Finally, a general relativistic 
description would be desirable, as the mass concentration 
and rotation rates of neutron stars clearly exceed the range 
for which a Newtonian treatment can be accurate. 
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